//$AVSTrans_TransFirst=0
//allow "else" instead of "," (for if()s)
//$AVSTrans_Replacement=else->,
//remove spaces
//$AVSTrans_Replacement= ->
//remove tabs
//$AVSTrans_Replacement=	->
//------------------------------------------------------------------
// these constants helped me to understand
C1S2=1/2;        // .5
C1S4=1/4;        // .25
C3S4=3/4;        // .75
C1S3=1/3;        // .3333333
C1S9=1/9;        // .1111111
C1S54=1/54;      // .0185185
CR3S2=sqrt(3)/2; // .8660254

x0r = 1000000;
x0i = 0;
x1r = 1000000;
x1i = 0;
x2r = 1000000;
x2i = 0;
x3r = 1000000;
x3i = 0;
if (equal(qa0, 0),
  if (equal(qb0, 0),
    if (equal(qc0, 0)
      // solving : qd0.x + qe0 = 0
      x0r = -qe0/qd0;
    else
      // solving : qc0.x^2 + qd0.x + qe0 = 0
      disc = sqr(qd0)-4*qc0*qe0;
      if (below(disc, 0),
        dr = 0;
	di = sqrt(-disc);
      ,
	dr = sqrt(disc);
        di = 0;
      );
      invsa = 1/(2*qc0);
      x0r = (-qd0+dr)*invsa;
      x0i = di*invsa;
      x1r = (-qd0-dr)*invsa;
      x1i = -di*invsa;
    );
  else // qb0 != 0
    // solving : qb0.x^3 + qc0.x^2 + qd0.x + qe0 = 0
    inv = 1/qb0;
    a2 = qc0*inv;
    a1 = qd0*inv;
    a0 = qe0*inv;
    q = (3*a1-sqr(a2))*C1S9;
    r = ((9*a1-2*sqr(a2))*a2-27*a0)*C1S54;
    disc = sqr(q)*q+sqr(r);
    /* x0r = 0; */
    /* x0i = 0; */
    /* x1r = 0; */
    /* x1i = 0; */
    /* x2r = 0; */
    /* x2i = 0; */
    if (below(disc, 0),
      theta = acos(r/sqrt(-sqr(q)*q));
      x0r = 2 * sqrt(-q) * cos(theta*C1S3) - a2*C1S3;
      x1r = 2 * sqrt(-q) * cos((theta+2*$PI)*C1S3) - a2*C1S3;
      x2r = 2 * sqrt(-q) * cos((theta+4*$PI)*C1S3) - a2*C1S3;
    else // disc >= 0    
      s = sign(r+sqrt(disc))*pow(abs(r+sqrt(disc)), C1S3);
      t = sign(r-sqrt(disc))*pow(abs(r-sqrt(disc)), C1S3);
      x0r = -a2*C1S3 + (s+t);
      x1r = -a2*C1S3 - (s+t)*C1S2;
      x2r = -a2*C1S3 - (s+t)*C1S2;
      x1i = CR3S2*(s-t);  /* abs(CR3S2*(s-t)) ?? */
      x2i = -x1i;
    );
  );
else // qa0 != 0
  // solving : qa0.x^4 + qb0.x^3 + qc0.x^2 + qd0.x + qe0 = 0
  inv = 1/qa0;
  qa = qb0*inv;
  qb = qc0*inv;
  qc = qd0*inv;
  qd = qe0*inv;
  a2 = -qb;
  a1 = qa*qc-4*qd;
  a0 = qd*(4*qb-sqr(qa))-sqr(qc);
  q = (3*a1-sqr(a2))*C1S9;
  r = ((9*a1-2*sqr(a2))*a2-27*a0)*C1S54;
  disc = sqr(q)*q+sqr(r);
  /* x0r = 0; */
  /* x0i = 0; */
  /* x1r = 0; */
  /* x1i = 0; */
  /* x2r = 0; */
  /* x2i = 0; */
  if (below(disc, 0),
    theta = acos(r/sqrt(-sqr(q)*q));
    /* x0r = 2*sqrt(-q)*cos(theta*C1S3)-a2*C1S3; */
    /* x1r = 2*sqrt(-q)*cos((theta+2*$PI)*C1S3)-a2*C1S3; */
    /* x2r = 2*sqrt(-q)*cos((theta+4*$PI)*C1S3)-a2*C1S3; */
    qz1 = 2 * sqrt(-q) * cos(theta*C1S3) - a2*C1S3; /* max(x0r,max(x1r,x2r)); */
  else // disc >= 0
    s = sign(r+sqrt(disc))*pow(abs(r+sqrt(disc)),C1S3);
    t = sign(r-sqrt(disc))*pow(abs(r-sqrt(disc)),C1S3);
    /* x0r = -a2*C1S3+s+t; */
    /* x1r = -a2*C1S3 -(s+t)*C1S2; */
    /* x2r = -a2*C1S3-(s+t)*C1S2; */
    /* x1i = abs(CR3S2*(s-t)); */
    /* x2i = -x1i; */
    qz1 = -a2*C1S3+s+t; /* x0r; */
  );
  qz2 = qz1+sqr(qa)*C1S4-qb;
  qz = sqrt(abs(qz2));
  if (equal(qz, 0),
    sq1 = C3S4*sqr(qa)-2*qb+2*sqrt(sqr(qz1)-4*qd);
    sq2 = C3S4*sqr(qa)-2*qb-2*sqrt(sqr(qz1)-4*qd);
  else // qz != 0
    sq1 = C3S4*sqr(qa)-qz2-2*qb+(4*qa*qb-8*qc-sqr(qa)*qa)/(4*qz);
    sq2 = C3S4*sqr(qa)-qz2-2*qb-(4*qa*qb-8*qc-sqr(qa)*qa)/(4*qz);
  );
  qv1 = sqrt(abs(sq1));
  qv2 = sqrt(abs(sq2));
  x0r = -C1S4*qa; /* x0r = 0; */
  x1r = -C1S4*qa; /* x1r = 0; */
  x2r = -C1S4*qa; /* x2r = 0; */
  x3r = -C1S4*qa; /* x3r = 0; */
  /* x0i = 0; */
  /* x1i = 0; */
  /* x2i = 0; */
  /* x3i = 0; */
  if (below(qz2, 0),
    x0i = C1S2*qz;
    x1i = C1S2*qz;
    x2i = -C1S2*qz;
    x3i = -C1S2*qz;
  else // qz2 >= 0
    x0r = x0r + C1S2*qz; /* x0r = C1S2*qz; */
    x1r = x1r + C1S2*qz; /* x1r = C1S2*qz; */
    x2r = x2r - C1S2*qz; /* x2r = -C1S2*qz; */
    x3r = x3r - C1S2*qz; /* x3r = -C1S2*qz; */
  );
  if (below(sq1, 0),
    x0i = x0i+C1S2*qv1;
    x1i = x1i-C1S2*qv1;
  else // sq1 >= 0
    x0r = x0r+C1S2*qv1;
    x1r = x1r-C1S2*qv1;
  );
  if (below(sq2, 0),
    x2i = x2i+C1S2*qv2;
    x3i = x3i-C1S2*qv2;
  else // sq2 >= 0
    x2r = x2r+C1S2*qv2;
    x3r = x3r-C1S2*qv2;
  );
);